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Abstract 

Understanding the earthquake (EQ) preparation process in terms of precursory electromagnetic (EM) emis- 
sions has been an evolving field of multi-disciplinary research. EM emissions in a wide frequency spectrum 
ranging from kHz to MHz are produced by opening cracks, which can be considered as precursors of gen- 
eral fracture. An important feature, observed on both laboratory and geophysical scale, is that the MHz 
radiation systematically precedes the kHz one. Yet, the link between an individual EM precursor and a 
distinctive stage of the EQ preparation comprises a crucial open question. A recently proposed two-stage 
model on preseismic EM activity suggests that the MHz EM emission is due to the fracture of the highly 
heterogeneous system that surrounds the fault. The finally emerged kHz EM emission is rooted in the final 
stage of EQ generation, namely, the fracture of entities sustaining the system. In this work we try to further 
penetrate and elucidate the link of the precursory kHz EM activity with the last stage of EQ generation 
building on two theoretical models for EQ dynamics. Firstly, the self-affine model states that an EQ is 
due to the slipping of two rough and rigid fractional Brownian profiles, one over the other, when there 
is an intersection between them. Secondly, the fragment- asperity model, rooted in a nonextensive Tsallis 
framework starting from first principles, consists of two rough profiles interacting via fragments filling the 
gap. In the latter approach, the mechanism of triggering EQ is established through the interaction of the 
irregularities of the fault planes and the fragments between them. This paper shows that these models of 
EQ dynamics can be linked with the detected kHz EM emission. In this framework of analysis of preseismic 
EM activity, we identify sufficient criteria that offer the possibility to discriminate whether a seismic shock 
is sourced in the fracture of fragments filling the gap between the rough profiles or in the fracture of "teeth" 
distributed across the fractional Brownian profiles that sustain the system. 

Keywords: Nonextensive Tsallis statistics, Preseismic Electromagnetic Emissions, Earthquake Dynamics, 
fractal Brownian motion model, self-affinity. Hurst exponent 



1. Introduction 

The use of basic principles of fracture mechanics for understanding the earthquake (EQ) preparation 
process is a challenging theme of multi-disciplinary research. Key fundamental questions in strength con- 
siderations of materials are: when does it fail? Are there signals that can warn of imminent failure? It has 
been shown that fracture induced electromagnetic (EM) physical fields allow real-time monitoring of damage 
evolution in materials during mechanical loading. EM emissions in a wide frequency spectrum ranging from 
kHz to MHz are produced by opening cracks, which can be considered as precursors of general fracture [1-9]. 
These precursors are detectable both at laboratory [1] and geophysical scale [10, 11]. An important feature, 
observed on both scales, is that the MHz radiation systematically precedes the kHz one [12-15]. Yet, the 
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link between an individual EM precursor and a distinctive stage of the EQ preparation comprises a crucial 
open question. 

The present paper draws on a recently proposed two-stage model on preseismic EM activity [14-19] 
suggesting that the MHz EM emission is due to the fracture of the highly heterogeneous system that 
surrounds the fault. This MHz EM activity can be attributed to phase transition of second order [15], while 
a Levy walk type mechanism can explain the observed critical state [16]. The finally emerged kHz EM 
emission, from approximately one week up to a few hours before the main shock occurrence, is rooted in the 
final stage of EQ generation, namely, the fracture of entities sustaining the system [14, 15, 17, 20, 21]. 

In this work we try to further penetrate and elucidate the link of the precursory kHz EM activity with the 
last stage of EQ generation building on two models for EQ dynamics. First, De Rubeis et al., (1996) [22] and 
Hallgass et al., (1997) [23] introduced a self-affine model for fault dynamics by means of the slipping of two 
rough and rigid Brownian profiles one over the other. In this self-affine asperity scheme, an individual EQ 
occurs when there is intersection between the two profiles. The second model concerns a fragment-asperity 
model for EQ dynamics which has been recently introduced by Solotongo-Costa and Posadas (SCP) [24]. 
This model is rooted in a nonextensive Tsallis framework starting from first principles, and draws on a 
scheme of two rough profiles interacting via fragments filling the gap. 

Herein, we argue that the aforementioned two models for EQ dynamics are identified in two qualitatively 
different epochs of kHz EM precursory activity: 

(i) the first epoch of the kHz EM precursor reflects the breakage of the fragments filling the gap between 
the two profiles. 

(ii) the second epoch is emitted during the fracture of large and strong entities ("teeth") which are dis- 
tributed along these two rough profiles of the fault. 

The structure of this paper is organized as follows: In section 2, we describe the aforementioned self- 
affine and fragment-asperity model for EQ dynamics. In Sec. 3 we set the context for the proposed model 
drawing on recent studies of kHz EM precursory activity in terms of complexity-organization. In Sec. 4 we 
analyze the kHz EM recordings by means of the Hurst's rescaled range analysis, distinguishing two emitted 
epochs that reflect on two different fracture regimes. In Sec. 5 we give the deflnition of the "electromagnetic 
earthquake" (EM-EQ) and we examine the nonextensive behaviour of the kHz EM activity focusing on the 
link between the two identified epochs associated with the last stage of EQ preparation process. Sec. 6 
draws on the footprints of fractal Brownian motion (fbm) profile and roughness of surface fracture related 
with the two epochs of kHz EM activity. In the latter part of that section we provide argumentation for 
the association of the second epoch with the self-affine model. Finally we summarize the key findings that 
support the proposed approach in view of further research considerations. 

2. Overview of two models for earthquake dynamics 

2.1. A self-ajfine asperity model for earthquake dynamics: experimental and theoretical evidence 

De Rubeis et al., (1996) [22] and Hallgass et al., (1997) [23] have worked on a model for regional fault 
dynamics by means of the slipping of two rough and rigid Brownian profiles one over the other. In this 
scheme, an individual EQ occurs when there is intersection between the two profiles. The authors assumed 
that the energy released is proportional to the extension of the overlap between the two asperities in contact. 
We briefly describe the rules of the model as defined by De Rubeis et al., (1996) [22]: 

(i) The initial condition is obtained by putting two rigid profiles in contact in the point where the height 
difference is minimal. 

(ii) The successive evolution is obtained by drifting a profile in a parallel way with respect to another. 

(iii) An intersection represents a single seismic event and starts with the collision of two asperities of the 
profiles. 
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This model exhibits a good interpretation of the scismicity generated in a large geographic area usually 
identified as "seismic region" , covering many geological faults, in a global sense. Ample experimental and 
theoretical evidence support the above mentioned scheme: 

(i) Kinematic or dynamic source inversions of EQs suggest that the final slip (or the stress drop) has a 
heterogeneous spatial distribution over the fault (see among others (author?) [25, 26, 27]). 

(ii) A study on Power spectrum analysis of the fault surface suggests that heterogeneities are observed 
over a large range of scale lengths [see Power et al., (1987), [28] in particular [Fig. 4]]. 

(iii) Investigators of the EQ dynamics have already pointed out that the fracture mechanics of the stressed 
crust of the earth forms self-similar fault patterns, with well-defined fractal dimensionalities [29-31]. 

(iv) Following the observations of the self-similarity in various length scales in the roughness of the fractured 
solid surfaces, Chakrabarti et al. (1999) [32] have proposed that the contact area distribution between 
two fractal surfaces follows a unique power law. 

(v) Huang and Turcottc, (1988) [33] pointed out that natural rock surfaces can be represented by fractional 
Brownian surfaces over a wide range. 

Interestingly, Hallgass et al. (1997) [23] have emphasized that "what is lacking is the description of 
what happened locally, i.e., as a consequence of a single event". In the following, we attempt to show that 
the statistics of regional seismicity could be merely a macroscopic reflection of the preparation process of a 
single EQ. 

2.2. A fragment- asperity model for earthquake dynamics coming from nonextensive statistical mechanics 

In seismology, the scaling relation between magnitude and the number of EQs is given by the Gutenberg- 
Richter (G-R) relationship: 

log N{> M) = a-bm. (1) 

where, A'^(> M) is the cumulative number of EQs with a magnitude greater than M occurring in a 
specified area and time and b and a are constants. It should be noted that the G-R empirical relation is not 
related with general physical principles. 

It is well-known that the Boltzmann-Gibbs statistical mechanics works best in dealing with systems 
composed of either independent subsystems or interacting via short-range forces, and whose subsystems can 
access all the available phase space. For systems exhibiting long-range correlations, memory, or fractal prop- 
erties, nonextensive statistical mechanics becomes the most appropriate mathematical framework [34. 35]. 
A central property of the EQ preparation process is the occurrence of coherent large-scale collective with a 
very rich structure, resulting from the repeated nonlinear interactions among its constituents. Consequently, 
the nonextensive statistical mechanics is an appropriate arena to investigate the last stage of EQ preparation 
process. 

As concerns the material between the fault planes, Herrmann and Roux (1990) [36] studied the phe- 
nomenon of fault slipping from a geometric viewpoint, offering an idealized representation of the fragmented 
core of a fault (gouge). They have presented the gouge as a self- medium formed by circular disk-shaped 
pieces which act like bearings filling the space between two planes [37]. In that direction, Sotolongo-Costa 
and Posadas (SCP) [24] have developed a model for EQ dynamics coming from a non-extensive Tsallis for- 
malism, starting from fundamental principles. In particular, the theoretical components of the SCP model 
read as follows: 

(i) The mechanism of relative displacement of fault plates is the main cause of EQs. 

(ii) The surfaces of the tectonic plates are irregular. The space between fault planes is filled with the 
residue of the breakage of the tectonic plates, from where the faults have originated. 
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(iii) The fragments are very diverse and have irregular shapes. The motion of the fault planes can be 
hindered not only by the overlapping of two irregularities/asperities (teeth) of the profiles, but also 
by the eventual relative position of several fragments. Thus, the mechanism of triggering EQs is 
established through the combination of the irregularities of the fault planes on one hand and the 
fragments between them on the other hand. 

(iv) In this regard, the authors studied the influence of the size distribution of fragments on the energy 
distribution of EQs. 

We note that the fragments size distribution function comes from a nonextensive Tsallis formulation, 
starting from first principles, i.e., a nonextensive formulation of the maximum entropy principle. Englaman 
et al. [38] showed that the standard Boltzmann-Gibbs formalism, although useful, cannot account for 
an important feature of the fragmentation process, i.e., the presence of scaling in the size distribution of 
fragments, which is one of the main ingredients of the SCP approach. Interestingly, the latter nonextensive 
approach leads to a G-R type law for the magnitude distribution of EQs (see Eq (8) in Ref. [24]). 

The fragment-asperity (SCP) model has been recently revised by Silva et al. [39] with two crucial 
updates. They use a different definition of the mean values in the context of Tsallis nonextensive statistics 
as proposed by Abe and Bagci [40]. Moreover, Silva et al. have introduced a new scale law, e (x r^, 
between the released energy e and the size r of fragments. The new scale proposed by Silva et al. [39] is 
in full agreement with the standard theory of rupture, namely, the well-known seismic moment scaling with 
rupture length (sec Ref. [41] for details). Finally, their approach leads to the following G-R type law for the 
magnitude distribution of EQs: 



log [N (> M)] = log iV + 

where, is the total number of EQs, N{> M) the number of EQs with magnitude larger than M, and 
M « log(e). The parameter a is the constant of proportionality between the EQ energy, e, and the size of 
fragment. The authors mentioned that this is not a trivial result, and incorporates the characteristics of 
nonextensivity into the distribution of EQs by magnitude. 

The entropic index q characterizes the degree of non-extensivity of the system reflected in the following 
pseudo-additivity rule: 



SgiA + B)= SM) + S<i{B) + (1 - q)S^{A)S^{B) (3) 

where, 

(4) 

is the Tsallis entropy, pi are the probabilities associated with the microscopic configurations, W is their 
total number, q is the nonextensive entropic parameter, and k is Boltzmann's constant. 

Importantly, from most geological analysis performed so far, values oi q ^ 1.6 — 1.8 seem to be universal, 
in the sense that different data sets from different regions of the globe indicate a value for the nonextensive 
parameter lying in this interval [24, 39, 41-43]. 



3. Linking earthquake dynamics and kHz electromagnetic activity: setting the context for the 
proposed approach 

Our analysis focuses on a well documented [14-17] kHz EM precursor (depicted in Fig. 1 with green 
and red color) associated with the Athens EQ that occurred on September 7, 1999 at 11:56 (GMT) with 
magnitude Mw = 5.9. In this paper, the preseismic magnetic field recorded by the 10 kHz North-South 
(NS) sensor is studied. 
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Drawing from recent studies, multidisciplinary analysis applied on recordings from the 10 kHz East- West 
(EW) sensor, has shown that the black part in Fig. 1 refers to the EM background activity (noise) which 
has been characterized by a low organization (or high complexity) [44, 45]. The first part of the precursor 
(epoch 1) has been characterized by the appearance of a population of EM events sparsely distributed in 
time with noteworthy higher order of organization in comparison to that of the noise. The strong bursts 
A and B in epoch 2 are characterized by a population of EM events of significantly higher organization in 
comparison to that of epoch 1 which are densely distributed in time. 



10 kHz NS : 28-Aug-1999 00:00:00 to 08-Sep-1999 00:00:00 




Time (days) 



Figure 1: View of the preseismic EM emission (green and red) recorded by the 10 kHz NS magnetic field 
sensor. The vertical arrow indicates the time of the Athens EQ occurrence. 



The same results are further verified here by the analysis of the recordings of the 10 kHz NS sensor. 
Characteristically in Fig. 2 we present the degree of organization between the two distinctive epochs namely 
1 and 2, by means of a robust grammar-based complexity/information technique, the T-Entropy. The 
analysis has been applied in terms of symbolic dynamics by using sequential successive windows of 1024 
samples each. We recall that T-entropy is based on the intellectual economy one makes when rewriting a 
string according to some rule. The estimation of T-Entropy has been presented in (author?) [21, 44]. 




Figure 2: Temporal evolution of T-Entropy using sequential blocks of 1024 samples each. 



The aforementioned results intuitively lead to the suspicion that the underlying fracture mechanisms 
are different in epoch 1 and epoch 2, characterized by a negative and positive-feedback mechanism, corre- 
spondingly. Against the theoretical background set in the previous section, herein we describe our proposed 
approach that aims to prove that the two models for EQ dynamics (described in Sec. 2) are identified in 
two qualitatively different epochs of the kHz EM precursory activity: 
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(i) The first epoch (Fig. 1, green color) of tire kHz EM precursor may reflect the breakage of the fragments 
filling the gap between the two profiles that sustain the fault (Fig. 3a). In this scheme, EM fluctuations 
are emitted during the fracture of fragments. 

(ii) The second epoch (Fig. 1, red color) is possibly emitted during the fracture of large and strong entities 
( "teeth" ) which are distributed along two rough and rigid fractal Brownian profiles one over the other 
(Fig. 3b). In this scheme, EM fluctuations are emitted during the fracture of strong and large "teeth". 



Figure 3: (a) An illustration of the fragment-asperity model. EM fluctuations are emitted during the fracture 
of fragments, (b) An illustration of the self-affine model. EM fluctuations are emitted during the fracture 
of strong and large "teeth" . 

In the next section we intend to verify that the two precursory epochs are rooted in different fracture 
regimes. For this purpose we first use the "Rescaled Range Analysis (R/S)" [46]. 

4. Analysis of presismic kHz EM emission in terms of the Hurst exponent 

After the works of Hurst et al. [46], Mandelbrot and Wallis [47] and Feder [48], Hurst's rescaled analysis 
has been used as a method to detect correlations in time series. There are two factors used in the R/S 
analysis: first, the range R, which is the difference between the maximum and minimum amounts of accu- 
mulated departure of the time series from the mean over a time span r, and second, the standard deviation 
S calculated over the time span r. The so-called rescaled range is exactly the ratio of R and S. From a 
variety of time series of natural phenomena, it is concluded that the ratio R/S is very well described by the 
following empirical relation: 



The correlation of the past and the future in the observational time series can be described by the H. 
For H = 0.5, there is an independent random process, with no correlations among samples. For H > 0.5, 
the sequence is characterized by a persistent behaviour, which means that the increasing or decreasing trend 
may followed by the same positive sign. For H < 0.5, the sequence is characterized by the anti-persistent 
behaviour, which means that an increasing or decreasing trend is more likely to be reversed implying a 
negative feedback mechanism. 

Focusing on the first epoch, we study the temporal evolution of the H exponent by applying a sequence 
of successive fixed time windows of 1024 samples each. The results obtained from Fig. 4, reveal that the 
first epoch is really characterized by anti-persistent behaviour {0 < H < 0.5) with a mean H = 0.38. 

Focusing on the second epoch, the same windowing method is used as that of epoch 1. Figs 5a and 5b 
depict the temporal evolution of H exponent for the two strong EM bursts contained in epoch 2. Results 
obtained from these figures reveal that both the two EM bursts arc characterized by a persistent (0.5 < 
H < 1) behaviour, with a mean H = 0.69 and H — 0.71 respectively. 




R/S = (t/2) 
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Figure 4: Temporal evolution of Hurst exponent (H) for epoch 1, deriving from fixed sequential windows 
of 1024 samples each. The green line is the threshold of the transition from antipersistent to persistent 
behavior. The pie shows that 97% of the blocks calculated, are antipersistent. 
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Figure 5: Temporal evolution of Hurst exponent (H) for the two strong EM bursts contained in epoch 2. 



To further support the persistency of the two strong EM bursts we applied the R/S method to the whole 
part of each one. Figs 6a and 6b depict the slopes deriving from the linear regression fitting of R/S method 
with H — 0.69 ± 0.05 and H = 0.71 ± 0.05 for the first and second EM bursts, respectively. The derived H 
exponents along with those of the windowing method strongly reveal that both EM bursts are characterized 
by an onset of persistency. 

Summarizing, the results of R/S analysis strongly support the hypothesis that the two preseismic epochs 
are rooted in different regimes characterized by negative and positive feedback correspondingly. To further 
investigate the possible differentiation of the two epochs we turn to the nonextensive model as described in 
section 2.2. 



5. Nonextensive approach of the two distinct epochs of kHz EM activity 

The analysis presented in this section first deals with the examination of whether the nonextensive 
formula 2 also describes the sequence of "electromagnetic earthquakes" (EM-EQs) included in the recorded 
EM precursor. Should this case exist, we go on to further examine the variations of the non-extensive 
parameter q and energy density a for different thresholds of magnitudes of the detected sequence of the 
EM-EQs. 

The definition of electromagnetic earthquake: Herein, we regard as amplitude ^ of a candidate "fracto- 
electromagnetic fluctuation" the difference Afj^^^iti) = A{ti) — Anoise, where Andse is the background (noise) 
level of the EM time series. We consider that a sequence of k successively emerged "fracto-electromagnetic 
fluctuations" Afj^j^j{ti), i = l,...,fc represents the damage of a fragment or a "tooth" in the irregular 
surfaces of fault. We refer to this as an "EM-EQ". Since the squared amplitude of the fracto-EM emissions 
is proportional to their energy Seai, the magnitude M of the candidate EM-EQ is given by the relation: 
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M = logesM - log [Af,,, {t,)f) (6) 

In order to achieve an efficient and accurate fit for the estimation of the nonextensive parameters q and 
a, the Levenberg-Marquardt (LM) method [49-53] was applied. The LM method is used to solve nonlinear 
least squares problems when fitting a parameterized function to a set of given data points. Fitting is achieved 
by minimizing the sum of the squares of the errors between the data points and the function. 

Figs. 7a and 7b show that Eq. 2 along with the use of LM method provides an excellent fit to the kHz 
EM-EQs included in the epochs 1 and 2, respectively, incorporating the characteristics of nonextensivity 
statistics into the detected kHz EM precursor. Herein, N is the total number of the detected EM-EQs, 
G(> M) = N{M >)/N the normalized cumulative number of EM-EQs with magnitude larger than M, 
and a the constant of proportionality between the EM energy released and the size of fragment [24, 39]. 
The best-fit parameters for this analysis are given by g = 1.739 ± 0.001 for epoch 1 and q — 1.834 ± 0.001 
for epoch 2 respectively. These values evidently reveal that both epochs 1 and 2 are characterized by high 
non-extensivity. 

It is very interesting to observe the similarity in the q values associated with Eq. 2 for all the catalogs 
of EQs used [q ~ 1.6 — 1.8, see Sec. 5), as well as for all the precursory sequences of EM-EQs under study. 
The observed similarity in the q values also indicates that the activation of a single EQ (fault) could be 
considered as a reduced self-affine image of the whole regional seismicity. These results are consistent with 
the self-affine nature of fracture and faulting processes [54]. 

We clarify that the parameter q itself is not a measure of the complexity-organization of the system but 
only a measure for the degree of non-extensivity. The dynamic changes of the complexity of the system 
can be quantified by applying the Tsallis entropy equation {Sq) to the time variation of the signal with a 
given nonextensive parameter q [20]. In that case, lower Sq values characterize the sequential portions of 
the signal with lower complexity (or higher organization). 

Figs. 8a and 8b depict the temporal variation of the Tsalis entropy for the epochs 1 and 2 respectively. 
The analysis has been applied in terms of symbolic dynamics by using sequential successive windows of 1024 
samples each. For each epoch, we use the estimated parameters q as derived from Eq. 2: q = 1.739 for 
epoch 1 and q = 1.834 for epoch 2. A detailed analysis of Tsallis entropy is presented in [20]. The results 
of this statistical analysis are consistent with those obtained from the T-Entropy in section 3. As shown, 
epoch 2 is characterized by a high organized population of EM events which are densely distributed in time 
in contrast to epoch 1, which is characterized by a population of lower organization sparsely distributed in 
time. Both epochs (1,2) have higher organization in respect to that of the backgroimd noise. 
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We recall that the scope of this analysis is to examine whether the two epochs of the precursory activity 
are really rooted in different regimes. The latter results support this hypothesis. In the following, we focus 
on the variation of the nonextensive parameter q and the volumetric energy density a for different cutoffs of 
magnitudes, Mc, of the detected EM-EQs, included in the two epochs of the emerged EM precursor, as the 
EQ is approaching. For statistically valid results, we used different Mc with an increment step of Mc = 0.1 
and a minimum event number of 50 events. For each step we fitted Eq. 2 to the available data points using 
the aforementioned LM method. In Figs. 9a and 9b the black curves depict the variation of the g-parameter 
for different M^- 

The observed high g-values for different threshold magnitudes reveal that the nonextensivity in the 
underlying fracture mechanism remains high in both epochs. This evidence is consistent with the hypothesis 
that the recorded EM precursor is associated with the final stage of EQ generation. Additionally, the 
prospective decrement of the nonextensive parameter q as the magnitude threshold increases is explained 
by the fact that the larger the magnitude threshold the larger the number of EM-EQs/fractures is omitted. 
Indeed, the large number of the small fractures along with the corresponding redistribution of stresses 
significantly contributes to the increment of the correlation length during the fracture process [55, and 
references therein]. The characteristic value that governs the overall system is the one that corresponds to 
the smaller magnitude threshold for each period. It should also be noted that although the nonextensive 
parameter q decreases at higher magnitude thresholds it still remains high. This evidence further verifies 
that strong correlations have developed within the system. 

We recall that a is the coefficient of proportionality between fragment size and released energy in the 
fragment-asperity model [39]. The blue curves in figs. 9a and 9b show that the energy density a characteristic 
value, increases at higher threshold values. This feature is consistent with the hypothesis that larger EM- 
EQs are rooted in larger and stronger entities. Special attention should be given to the following points of 
differentiation between the two highly no n- nonextensive epochs 1 and 2: 

(i) The maximum magnitude of EM-EQs in epoch 1 is M„iax = 8.19 while in epoch 2 it is M„iax ~ 10.36. 

(ii) EM-EQs in the abruptly emerged epoch 2 reach significantly higher values of energy density a in the 
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Figure 8: Temporal evolution of Tsallis entropy for epochs 1 and 2 respectively. The blue line depicts the 
entropy level that refers to the background EM activity. 
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Figure 9: Variation of nonextensive parameter q and the volumetric energy density a, for different thresholds 
of magnitudes of the detected EM-EQs for epoch 1 and epoch 2 respectively. 



order of 10^^ while in epoch 1 reach values in the order of 10^^. 
Summarizing, up to this stage of the analysis we have shown that: 

(i) The degree of the organization of the two strong EM bursts included in epoch 2 is significantly higher 
than in epoch 1. 

(ii) The abrupt emerged two kHz EM bursts contained in epoch 2 are characterized by persistency in 
contrast to epoch 1 which is characterized by antipersistency. 

(iii) Both two epochs arc characterized by strong noncxtensivity. 

(iv) The second epoch is rooted in the fracture of larger and stronger entities in comparison to the first 
epoch. 

The above mentioned characteristics imply that the two epochs of prescismic kHz EM activity under study 
refer to different models of EQ dynamics. In the following we attempt to establish the hypothesis that the 
first and second epochs arc consistent with the nonextensive fragment and sclf-affinc model, correspondingly. 
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6. Footprints of persistent-fBm profile of fracture surfaces in the second epoch of kHz EM 
activity 

Despite the complexity of the Earth's crust there are several universally holding scaling relations [21, and 
references therein] . Such universal structural patterns of fracture and faulting process should be included 
into an associated EM precursor. Notice that Maslov et al. (1994) [56] have formally established the 
relationship between spatial fractal behavior and long-range temporal correlations for a broad range of 
critical phenomena. They showed that both the temporal and spatial activity can be described as different 
cuts in the same underlying fractal. 

From the early work of Mandelbrot (1982) [57], much effort has been put to statistically characterise the 
resulting fractal surfaces in fracture processes: 

1. Fracture surfaces were found to be self-afhne following the fractional Brownian motion (fBm) model 
over a wide range of length scales. 

2. The spatial roughness of fracture surfaces has been interpreted as a universal indicator of surface 
fracture, weakly dependent on the nature of the material and on the failure mode. 

Therefore, a fracture surface follows the persistent fBm-model, and consequently, an associated EM 
precursor should behave as a persistent fBm temporal fractal. 

Herein, we note that the characteristic of the persistent fBm model does not refer to a population of 
residual-fragments coming from previous fractures. Consequently, if our hypothesis is true, this characteristic 
should be contained only in epoch 2 due to the interference of the large and strong teeth. 

6.1. Footprints of universal roughness of fracture surfaces 

The Hurst Exponent H specifies the strength of the irregularity ( "roughness" ) of the surface topography 
[58] and it has also been interpreted as a universal indicator of fracture surfaces [59-62]. The height- height 

1 /2 

correlation function of a surface A/i(r) =< [h{r -\- Ar) — h(r)\r > computed along a given direction has 
been found to scale as Ah ~ (Ar) , where H refers to the Hurst exponent. The Hurst exponent H ^ 0.7 has 
been interpreted as a universal indicator of surface fracture, weakly dependent on the nature of the material 
and the failure mode [59-63]. Importantly, Renard et al., (2006) [64] measured the surface roughness of a 
recently exhumed strike-slip fault plane by three independent 3D portable laser scanners. Their statistical 
scaling analyses revealed that the striated fault surface exhibits self-affine scaling invariance that can be 
described by a scaling roughness exponent. Hi = 0.7 in the direction of slip. In Section 4 we showed that 
the "roughness" of the profile of the two strong kHz EM bursts, as it is represented by the Hurst exponent, 
is distributed around the value of 0.7. Thus, the universal spatial roughness of fracture surfaces nicely 
coincides with the roughness of the temporal profile of the recorded two strong preseismic kHz EM bursts 
included in the second epoch and does not coincide with the temporal profile of the preseismic time series 
of the first epoch. Herein, we recall the mean Hurst exponent of epoch 1 is H = 0.38. 

6.2. Footprints of fractional- Brownian-motion model of fracture surfaces 

If a pre-seismic EM time series behaves as a persistent temporal fractal, then, a power-law of the form 
S{f) oc f~^ is obeyed, with S{f) the power spectral density (PSD) and / the frequency. In a log 5(7) — log/ 
representation the power spectrum is a line with slope /3. The spectral scaling exponent /3 is a measure of 
the strength of time correlations, while the quality of the fit of a time series to the power-law is represented 
by the linear correlation coefficient r. 

The "wavelet spectrum" is used in order to provide an unbiased and consistent estimation of the true 
power spectrum of the time-series with the 'Morlet' wavelet as a mother function. A fixed 1024 moving 
window with no overlapping sequence was used for the bursts, 2(A) and 2(B), see Fig. 1 respectively. For 
each window the local parameters /3 and r were derived. In Figs. 10a and 10b, we present the distribution 
of the /3- exponents with r > 0.95. 
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Figure 10: The distribution of the (3 - exponents with r > 0.95 as apphed to epoch 2(A) and 2(B) respectively 

Two classes of signal have been widely used to model stochastic fractal time series: fractional Gaussian 
noise (fGn) and fractional Brownian motion (fBm) [65] . For the case of the fCn model the scaling exponent 
/3 lies between —1 < /3 < 1, while the regime of fBm is indicated by /3 values from 1 to 3. The distribution of 
/3 - exponents reveal that the profile of the two strong EM bursts follows the fBm-model. To further support 
this finding, Figs. 11a and lib show the linear regression fitting for the estimation of the /3 exponents of 
the overall two strong EM bursts with values: /3 = 2.28 ± 0.06 and /3 = 2.43 ± 0.09 for the epochs 2(A) 
and 2(B). This finding verifies that the profile of the second epoch of the candidate kHz EM precursor is 
qualitatively analogous to the fBm-model (1 < /3 < 3), and could be originated during the slipping of two 
rough and rigid Brownian profiles [22] that follow the fBm-model. 



Waveiet Spectra Anaiysis, w=1024 
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Figure 11: Overall linear regression fitting for the estimation of the /? exponent of the power spectra density 
function: S{f) oc /"^ 



We recall that the /3 exponent is related to the Hurst exponent, H, by the formula /3 = 2H + 1 with 
< if < 1 (1 < /3 < 3) for the fBm model [65]. Fig. 11 shows the estimated values of Hurst exponent, 
H = 0.64 and H = 0.72, calculated by applying the (3 - exponent to the formula: /? = 2H + 1. This finding 
further verifies that two strong EM bursts follow the persistent fBm-model. 

Fig. 12 depicts the distribution of iJ-exponents deriving from the formula (3 = 2H -I- 1 in relation to the 
distribution of iJ-cxponents as they have been estimated in terms of R/S analysis, for the epochs 2(A), 2(B) 
and 1, respectively. We observe that the two distributions (Figs. 12a and 12b) referred to the two strong 
EM bursts contained in epoch 2 concur, supporting the persistcnt-fBm temporal profile of the second epoch. 
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On the contrary the iJ-exponcnt distributions of epoch 1 (Fig. 12c), verify the antipersistcnt profile which 
is not consistent with the universal indicator of fracture surfaces {H ~ 0.7). 
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Figure 12: The distribution of iy-exponents deriving from the formula /3 = 2H + 1 in relation to the 
distribution of Hurst /i-expoucnts as they have been estimated in terms of R/S analysis, respectively from 
epoch 2(A), 2(B) and epoch 1. 



6.3. Additional arguments supporting the association of the second epoch with the self-affine asperity model 

In this section we refer to arguments which further support the association of the second epoch with the 
self-affine asperity model introduced by De Rubeis et al. [22] and Hallgass et al. [23] . 

According to the self-afhne asperity model, the distribution of areas of the asperities broken A, follows 
a power law P{A) ^ A~^ , with an exponent 5, which could be related to the Hurst exponent (0 < < 1) 
that controls the roughness of the fault. The former relation is obtained by supposing that the area of the 
broken asperities scales with its linear extension / as Aasp ^ V-^^ . The H-cxponent associated with the 
second epoch of the kHz EM activity is distributed around the value 0.7. It is reasonable to assume that the 
broken teeth scaled with its linear dimension / as Aasp ^ Z^ ''^. Importantly, numerical studies performed 
by de Arcangelis et al. (1989) [66] indicate that the number of bonds that break, scales during the whole 
process of fracture as 1.7. This consistency further supports the association of the second epoch with the 
self-affine asperity model. 

The self-affine model also reproduces the Gutenberg-Richter law suggesting that a seismic event releases 
energy m the interval [E, E+dE] with a probability P{E)dE, P{E) - E'^ , where B = a+1 and a = l-H/2 
with a G [1/2, 1]. In the present case, the Hurst-exponent H ~ 0.7 leads to a ~ 0.65. Thus, the fracture of 
asperities released EM energies following the distribution P{E) ~ E~^ , where B ^ 1.65. This value is in 
harmony with geophysical data. Indeed, the distribution of energies released at any EQ is described by the 
power-law, P{E) - E'^, where B - 1.4 - 1.6 [67]. 

7. Discussion &: Conclusions 

EM emissions in a wide frequency spectrum ranging from kHz to MHz are produced by opening cracks, 
which can be considered as precursors of general fracture. According to a recently proposed two-stage 
model on preseismic EM activity, the MHz EM emission is thought to be due to the fracture of the highly 
heterogeneous system that surrounds the fault. The finally emerged kHz EM emission is rooted in the final 
stage of EQ generation, namely, the fracture of entities sustaining the system. In this work we have further 
examined and elucidated the link of the precursory kHz EM activity with the last stage of EQ generation 
building on two theoretical models for EQ dynamics. The first sclf-affinc asperity model states that the 
EQ is due to the slipping of two rough and rigid Brownian profiles one over the other. In this scheme, 
an individual EQ occurs when there is an intersection between the two fractional Brownian profiles. The 
second model, which is rooted in a nonextensive Tsallis framework starting from first principles, consists of 
two rough profiles interacting via fragments filling the gap. In this nonextensive approach, the mechanism of 
triggering EQ is established through the interaction of the irregularities of the fault planes and the fragments 
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between them. This paper has attempted to show that two models of EQ dynamics supplement each other, 
in a sense, and both are mirrored in the detected two qualitative different epochs of the preseismic kHz 
EM emission associated with the Athens EQ. We have argued that the initially emerged precursory activity 
(epoch 1) follows the nonextensive model and is due to the fracture of fragments filling the gap between the 
two rough planes of the activated fault. The finally emerged precursory activity containing the two abruptly 
emerged strong impulsive EM bursts, follows the self-affine asperity model. 

More precisely, our results illustrate that both epochs are characterized by strong nonextensivity. This 
finding supports the association of two epochs with the final stage of EQ generation. However, significant 
differentiations between the two distinct epochs of the recorded kHz EM activity, verify the link with the 
two models of EQ dynamics: 

(i) The degree of the organization of the two strong EM bursts included in epoch 2 is significantly higher 
than in the initial epoch 1. 

(ii) Importantly, the abrupt emerged two kHz EM bursts in the tail of the preseismic EM activity (epoch 
2) are characterized by persistency in contrast to epoch 1 which is characterized by antipersistency. 

Recently, Sammis and Sornette (2002) [68] presented the most important mechanisms for such positive 
feedback. Laboratory experiments by means of acoustic and EM emission show that the main rupture 
occurs after the appearance of persistent behaviour ([69-71]. Recently Lei et al. (2000, 2004) [72, 73] 
have studied how individual and coupled asperities fracture as well as the role of asperities in fault 
nucleation and as potential precursors prior to dynamic rupture. The self-excitation strength, which 
expresses the strength of the effect of excitation associated with the preceding event on succeeding 
events, or equivalently, the degree of positive feedback in the dynamics, reaches the maximum value of 
1 during the nucleation. 

(iii) The second epoch is rooted in the fracture of larger and stronger entities in comparison to the first 
epoch. The maximum magnitude of EM-EQs in epoch 1 is M^ax ~ 8.19 while in epoch 2 it is 

Mmax = 10.36. 

(iv) EM-EQs in the abruptly emerged epoch 2 reach significantly higher values of energy density a of the 
order of 10^ while in epoch 1 energy reaches values up to 10^^. We further emphasise that the breakage 
of fragments filling the gap is easier than the breakage of teeth distributed along the rough fault planes 
[24]. 

(v) The kHz EM precursor of the second epoch includes crucial universal feature of fracture of surfaces; 
namely, it follows the persistent fBm-model with a roughness consistent with the universal roughness 
of fracture surfaces {H ~ 0.7). This universal footprint is not mirrored in the first epoch. 

The key differentiations mentioned above support the following hypothesis: the first epoch refers to the 
fracture of fragments intervening between the two anomalous surfaces of the fault that contributes to the 
hindering of their relative motion. Once the fracture of one fragment has occurred, there is a reformation of 
fragments followed by a redistribution of stresses. This process practically results to a relative displacement 
of the fault planes (fault slip). The next EM-EQ will emerge when a new fragment breaks under the 
impact of the increased tensions. This process is consistent with the antipersistent character of the first 
epoch. As the fragments are broken, the two rough planes of the fault approaches each other. The abruptly 
emerged second epoch of large EM-EQs is the reflection of collision and breakup of large and strong teeth of 
the irregular surfaces. Its persistent behaviour and high organization along with the corresponding higher 
magnitudes M and energy density a values supports this hypothesis. 

A challenging issue in the field of EQ dynamics is whether the generation of an EQ can be adequately 
explained by the self-afHne model described in section 2.1 or the nonextensive fragment-asperity model as 
described in section 2.2. Confronted with such a question, we may not seek answers merely in the statistics 
of EQs. Both two models of EQ dynamics lead to a power-law distribution of magnitudes. On the one hand, 
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the sclf-asperity model leads to the traditional power-law distribution of magnitudes, known as Gutenberg- 
Richter law. On the other hand, it has been shown [74] that the fragment-asperity model as described by 
Eq. 2 also leads to a Gutenberg- Riehter- type - law for large EQs by the following equation [74]: 



The fact that both models find quantitative expression through the latter equation is not an unexpected 
result. Indeed, some generic behaviour seems to be shared by fragmenting systems whatever their size, 
material, or typical interaction energy. For instance, as observed in a large variety of experiments [75, 76] 
and natural phenomena [76-79], the fragment size distribution frequently exhibits a power law behavior, the 
origin of which is still unknown. Several examples of power-law fragmentation are given in [80] . Considering 
that the released energy e is proportional to fragment size [39] , it is reasonable that the magnitudes of EQs 
which are rooted in the fracture of the population of fragments filling the space between fault planes also 
follow a power-law distribution. 

In the research prospect of discriminating whether a seismic shock is sourced in the fracture of fragments 
filling the gap between the rough profiles or in the fracture of teeth distributed across the fractional Brownian 
profiles, it would be more appropriate to focus on the generation of a single EQ and not on the statistics 
of a population of different EQs. The reported significant differentiations of the two epochs imply that 
the fracture induced kHz pre-seismic EM emissions seem to offer such a possibility. Our proposed approach 
cannot claim any universal applicability or objective truth. The experimental results presented in this paper 
justify the link between nonextensive fragment-asperity and self-affine asperity models for EQ dynamics 
with preseismic kHz EM activity. However, future research needs to be undertaken in order to test the 
applicability of such an approach in different data-sets and further verify the key features identifying with 
epoch 1 and epoch 2 respectively, as proposed herein. 
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